      SUBROUTINE VEGRAD(jmon,cos_z, tpbl, snow, asnow, fwet, albsod,& 
      & solad, solai, laisun, laisha, parsun, parsha)
!INPUT

      use INMSOIL_NUMPARAMS, only : &
     & ML, MS
      USE CBALANCE_CONST
      !INCLUDE 'soil.inc'
!      PARAMETER (nv2=7)
     ! USE PTRMODPHY, ONLY : NV2 
      INTEGER, PARAMETER :: twoband=2

      integer jmon !month number
	  REAL (KIND=8)  cos_z !cosine of direct zenith angle
	  REAL (KIND=8)  tpbl            !air temperature (kelvin)
	  REAL (KIND=8)  snow              !snow depth, m
	  REAL (KIND=8)  asnow    !snow albedo (direct)
	  REAL (KIND=8)  fwet(nv2) !wet franction of canopy (including snow)
!	  REAL (KIND=8)  wll !volumetric water content in upper soil layer
!	  REAL (KIND=8)  color_ich(9)	!fractions covered by different soil colors
	  REAL (KIND=8)  solad(twoband)       !incodmin1g direct solar radiation (w/m**2) 
      REAL (KIND=8)  solai(twoband)       !incodmin1g diffuse solar radiation (w/m**2) 
	  REAL (KIND=8)  albsod(twoband)      !ground albedo (direct)
      REAL (KIND=8)  albsoi(twoband)      !ground albedo (diffuse)



!OUTPUT
     
      REAL (KIND=8)  laisun(nv2)        !sunlit leaf area index, one-sided
      REAL (KIND=8)  laisha(nv2)        !shaded leaf area index, one-sided
      REAL (KIND=8)  parsun(nv2)        !average absorbed par for sunlit leaves (w/m**2)
      REAL (KIND=8)  parsha(nv2)        !average absorbed par for shaded leaves (w/m**2)

!LOCAL VARIABLES


	  REAL (KIND=8)  fsun(nv2)          !sunlit fraction of canopy
      REAL (KIND=8)  vai(nv2)           !leaf + stem area, one-sided
	  REAL (KIND=8)  elai(nv2)
	  REAL (KIND=8)  esai(nv2)
      REAL (KIND=8)  tv(nv2)            !vegetation temperature (kelvin)

      REAL (KIND=8)  fabd(twoband,nv2)    !flux abs by veg (per unit incodmin1g direct flux) 
      REAL (KIND=8)  fabi(twoband,nv2)    !flux abs by veg (per unit incodmin1g diffuse flux)
      REAL (KIND=8)  ftdd(twoband,nv2)    !down dir flux below veg (per incodmin1g dir flux) 
      REAL (KIND=8)  ftid(twoband,nv2)    !down dif flux below veg (per incodmin1g dir flux) 
      REAL (KIND=8)  ftii(twoband,nv2)    !down dif flux below veg (per incodmin1g dif flux) 
	  REAL (KIND=8) ftdi(twoband,nv2) !down direct flux below veg per unit dif flux = 0

      REAL (KIND=8)  albd(twoband,nv2)    !overall surface albedo (direct)
      REAL (KIND=8)  albi(twoband,nv2)    !overall surface albedo (diffuse)

	  REAL (KIND=8)  igs(nv2)
!	  REAL (KIND=8)  inc, albedo_dry, albedo_sat
!	  REAL (KIND=8)  albsat(9,2), albdry(9,2)
	  REAL (KIND=8)  rho(2,nv2),tau(2,nv2)
	  REAL (KIND=8)  wl, ws, ol, fb, ext, wl2
	  REAL (KIND=8)  prevent
	  REAL (KIND=8)  gdir(nv2)
      

  
      prevent = 0.000001


      rhol(1,1:nv2) = (/0.10,0.10,0.07,&
     &                          0.11,0.10,0.11,0.11/)
                              
     rhol(2,1:nv2) = (/0.45,0.45,0.35,&
     &                          0.58,0.45,0.58,&
     &                          0.58/) !,0.00/

! stem reflectance: 1=vis, 2=nir
      rhos(1,1:nv2) = (/0.16,0.16,0.16,&
     &                          0.36,0.16,0.36,&
     &                          0.36/) !,0.00/
      rhos(2,1:nv2) = (/0.39,0.39,0.39,&
     &                          0.58,0.39,0.58,&
     &                          0.58/) !,0.00/

! leaf transmittance: 1=vis, 2=nir
      taul(1,1:nv2) = (/0.05,0.05,0.05,&
     &                          0.07,0.05,0.07,&
     &                          0.07/) !,0.00/
      taul(2,1:nv2) = (/0.25,0.25,0.10,&
     &                          0.25,0.25,0.25,&
     &                          0.25/) !,0.00/

! stem transmittance: 1=vis, 2=nir
      taus(1,1:nv2) = (/0.001,0.001,0.001,&
     &                          0.220,0.001,&
     &                          0.220,0.220/) !,0.000/
      taus(2,1:nv2) = (/0.001,0.001,0.001,&
     &                          0.380,0.001,&
     &                          0.380,0.380/) !,0.000/
     open(36, file='source/include/tai.txt')
     open(37, file='source/include/gai.txt')
     read(36,*)((tai(j,i),i=1,12),j=1,nv2)
     read(37,*)((gai(j,i),i=1,12),j=1,nv2)
     close (36)
     close(37)

         do k = 1, nv2

            ol = dmin1( dmax1(snow-hbot(k),0.), htop(k)-hbot(k))
            fb = 1. - ol / dmax1(1.e-06, htop(k)-hbot(k))
            elai(k) = gai(k,jmon)*fb
            esai(k) = (tai(k,jmon)-gai(k,jmon))*fb
            if (elai(k) .lt. 0.05) elai(k) = 0.
            if (esai(k) .lt. 0.05) esai(k) = 0.
         end do
! -----------------------------------------------------------------
! weight reflectance/transmittance by lai and sai
      
         do i = 1, twoband
				do j = 1, nv2
					vai(j) = elai(j) + esai(j)
					wl = elai(j) / dmax1(vai(j),prevent)
					ws = esai(j) / dmax1(vai(j),prevent)
					rho(i,j)=dmax1(rhol(i,j)*wl+rhos(i,j)*ws, prevent)
					tau(i,j)=dmax1(taul(i,j)*wl+taus(i,j)*ws, prevent)
					ftdi(i,j)=0.0
				end do
			end do


          	 do i=1, twoband
			 albsoi(i)=albsod(i)
!---------------------------------------------------------------------
! if  soil is covered by snow. Added by AYu
!---------------------------------------------------------------------
              if (snow.NE.0.0)then   
				 albsod(i) = asnow
				 albsoi(i) = asnow  
			  end if 
!------------------------------------------------------------------------			    
		   
            end do

!  vegetation: solar fluxes for unit incodmin1g direct (ic=0) and diffuse flux (ic=1)

			do i = 1, twoband

				ic = 0
				call albveg (tpbl, fabd   ,albd   ,ftdd   ,ftid   ,albsod ,& 
     &                         albsoi ,i      ,ic     ,cos_z  ,& 
     &                         vai    ,fwet   ,rho    ,& 
     &                         tau, gdir)
				ic = 1
				call albveg (tpbl, fabi   ,albi   ,ftdi   ,ftii   ,albsod ,& 
     &                         albsoi ,i      ,ic     ,cos_z  ,& 
     &                         vai    ,fwet   ,rho    ,& 
     &                         tau, gdir)
			end do

! sunlit fraction of canopy. set fsun = 0 if fsun < 0.01.

			do i = 1, nv2
                        ! write(*,*) gdir(i), cos_z, vai(i)
				ext = gdir(i)/dmax1(cos_z,prevent) *&
     &				sqrt(1.-rho(1,i)-tau(1,i))
				fsun(i) = (1.-exp(-ext*vai(i))) / &
     &				dmax1(ext*vai(i),prevent)
				ext = fsun(i)                                 !temporary fsun
				if (ext .lt. 0.01) then 
					wl2 = 0.                                   !temporary fsun
				else
					wl2 = ext                                  !temporary fsun
				end if
				fsun(i) = wl2
			end do
     
	 call surrad (fsun, elai, esai, vai, solad, solai, fabd, fabi, ftdd, &
	 & ftid, ftii, albsod, albsoi, laisun, laisha, parsun, parsha)


      return
      end

!--------------------------------------------------------------------------------

		          
      subroutine albveg (tpbl, fab    ,fre    ,ftd   ,fti    ,albsod ,&
     &                   albsoi ,ib     ,ic    ,cos_z  ,&
     &                   vai    ,fwet   ,rho   ,&
     &                   tau ,gdir )
					   


!------------------------ code history ---------------------------
!source file:       twostr.F
!purpose:           two-stream fluxes for canopy radiative transfer
!date last revised: March 1996 - lsm version 1
!author:            Gordon Bonan
!standardized:      J. Truesdale, Feb. 1996
!reviewed:          G. Bonan, Feb. 1996
!-----------------------------------------------------------------
 
!------------------------ notes ----------------------------------
!use two-stream approximation of Dickinson (1983) Adv Geophysics
!25:305-353 and Sellers (1985) Int J Remote Sensing 6:1335-1372
!to calculate fluxes absorbed by vegetation, reflected by vegetation,
!and transmitted through vegetation for unit incodmin1g direct or diffuse 
!flux given an underlying surface with known albedo.
!-----------------------------------------------------------------

!------------------------ input/output variables -----------------

!input

      use INMSOIL_NUMPARAMS, only : &
      & ML, MS
      USE CBALANCE_CONST
       !INCLUDE 'soil.inc'
!      PARAMETER (nv2=13)
!      USE PTRMODPHY, ONLY : NV2 

      PARAMETER (twoband=2)

      integer ib              !waveband number 
      integer ic              !0=unit incodmin1g direct; 1=unit incodmin1g diffuse

      REAL (KIND=8)  tpbl
      REAL (KIND=8)  cos_z              !cosine of direct zenith angle
      REAL (KIND=8)  vai(nv2)           !one-sided leaf+stem area index
      REAL (KIND=8)  fwet(nv2)          !fraction of lai, sai that is wetted

      REAL (KIND=8)  albsod(twoband)    !direct  albedo of underlying surface 
      REAL (KIND=8)  albsoi(twoband)    !diffuse albedo of underlying surface
      REAL (KIND=8)  rho(twoband,nv2)   !leaf+stem reflectance
      REAL (KIND=8)  tau(twoband,nv2)   !leaf+stem transmittance


!output
      REAL (KIND=8)  fab(twoband,nv2)  !flux abs by veg layer (per unit incodmin1g flux)
      REAL (KIND=8)  fre(twoband,nv2)  !flux refl above veg layer (per unit incodmin1g flux)
      REAL (KIND=8)  ftd(twoband,nv2)  !down dir flux below veg layer (per unit in flux)
      REAL (KIND=8)  fti(twoband,nv2)  !down dif flux below veg layer (per unit in flux)
      REAL (KIND=8)  gdir(nv2)         !relative projected leaf+stem area in solar direction
!-----------------------------------------------------------------

!------------------------ local variables ------------------------
      REAL (KIND=8)  omega(nv2)   !fraction of intercepted radiation that is scattered
      REAL (KIND=8)  omegal(nv2)  !omega for leaves
      REAL (KIND=8)  betai(nv2)   !upscatter parameter for diffuse radiation 
      REAL (KIND=8)  betail(nv2)  !betai for leaves
      REAL (KIND=8)  betad(nv2)   !upscatter parameter for direct beam radiation 
      REAL (KIND=8)  betadl(nv2)  !betad for leaves
      REAL (KIND=8)  ext(nv2)     !optical depth of direct beam per unit leaf area 
      REAL (KIND=8)  avmu(nv2)    !average diffuse optical depth
      REAL (KIND=8)  tv(nv2)            !vegetation temperature (kelvin)
 
      integer j,i       !array index 

      REAL (KIND=8)  cosz         !0.001 <= coszen <= 1.000
      REAL (KIND=8)  asu          !single scattering albedo
      REAL (KIND=8)  chil         ! -0.4 <= xl <= 0.6

      REAL (KIND=8)  tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9
      REAL (KIND=8)  p1,p2,p3,p4,s1,s2,u1,u2,u3
      REAL (KIND=8)  b,c,d,d1,d2,f,h,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
      REAL (KIND=8)  phi1,phi2,sigma
      REAL (KIND=8)  ftds,ftis,fres
	  REAL (KIND=8)  tfrz, tgr
	  REAL (KIND=8)  omegas(twoband)



! omega,betad,betai for snow 
      omegas(1) = 0.8
      omegas(2) = 0.4


!-----------------------------------------------------------------

!calculate two-stream parameters omega, betad, betai, avmu, gdir, ext.
!omega, betad, betai are adjusted for snow. values for omega*betad 
!and omega*betai are calculated and then divided by the new omega
!because the product omega*betai, omega*betad is used in solution. 
!also, the transmittances and reflectances (tau, rho) are linear 
!weights of leaf and stem values.

	tfrz = 273.15

	do i = 1, nv2

!		fwet(i) = 0.5
		cosz = dmax1(0.001, cos_z)
		chil = dmin1( dmax1(xl(i), -0.4), 0.6)
		if (abs(chil) .le. 0.01) chil = 0.01
		phi1 = 0.5 - 0.633*chil - 0.330*chil*chil
		phi2 = 0.877 * (1.-2.*phi1)
		gdir(i) = phi1 + phi2*cosz
		ext(i) = gdir(i)/cosz
		avmu(i) = ( 1. - phi1/phi2 * log((phi1+phi2)/phi1) ) / phi2
		omegal(i) = rho(ib,i) + tau(ib,i)
		tmp0 = gdir(i) + phi2*cosz
		tmp1 = phi1*cosz
		asu = 0.5*omegal(i)*gdir(i)/tmp0 * ( 1. - tmp1/tmp0 *&
     &         log((tmp1+tmp0)/tmp1) )
		betadl(i) = (1.+avmu(i)*ext(i))/(omegal(i)*avmu(i)*ext(i))*asu
		betail(i) = 0.5* ( rho(ib,i)+tau(ib,i) + (rho(ib,i)-tau(ib,i)) &
     &               * ((1.+chil)/2.)**2 ) / omegal(i)
	end do

!adjust omega, betad, and betai for intercepted snow

      do i = 1, nv2
	     tv(i)=tpbl
	     tgr = tv(i) + 3.   !AYu
         if (tgr .gt. tfrz) then                                !no snow
            tmp0 = omegal(i)           
            tmp1 = betadl(i) 
            tmp2 = betail(i)  
         else
            tmp0 = (1.-fwet(i))*omegal(i) + fwet(i)*omegas(ib)           
            tmp1 = ( (1.-fwet(i))*omegal(i)*betadl(i) +&   
     &              fwet(i)*omegas(ib)*betads ) / tmp0
            tmp2 = ( (1.-fwet(i))*omegal(i)*betail(i) +&  
     &              fwet(i)*omegas(ib)*betais ) / tmp0
         end if
         omega(i) = tmp0           
         betad(i) = tmp1 
         betai(i) = tmp2  
      end do

!absorbed, reflected, transmitted fluxes per unit incodmin1g radiation

      do i = 1, nv2
         b = 1. - omega(i) + omega(i)*betai(i)
         c = omega(i)*betai(i)
         tmp0 = avmu(i)*ext(i)
         d = tmp0 * omega(i)*betad(i)
         f = tmp0 * omega(i)*(1.-betad(i))
         tmp1 = b*b - c*c
         h = sqrt(tmp1) / avmu(i)
         sigma = tmp0*tmp0 - tmp1
         if(sigma.eq.0.) sigma = 0.0001
         p1 = b + avmu(i)*h
         p2 = b - avmu(i)*h
         p3 = b + tmp0
         p4 = b - tmp0
         s1 = exp(-h*vai(i))
         s2 = exp(-ext(i)*vai(i))
         if (ic .eq. 0) then
            u1 = b - c/albsod(ib)
            u2 = b - c*albsod(ib)
            u3 = f + c*albsod(ib)
         else
            u1 = b - c/albsoi(ib)
            u2 = b - c*albsoi(ib)
            u3 = f + c*albsoi(ib)
         end if
         tmp2 = u1 - avmu(i)*h
         tmp3 = u1 + avmu(i)*h
         d1 = p1*tmp2/s1 - p2*tmp3*s1
         tmp4 = u2 + avmu(i)*h
         tmp5 = u2 - avmu(i)*h
         d2 = tmp4/s1 - tmp5*s1
         h1 = -d*p4 - c*f
         tmp6 = d - h1*p3/sigma
         tmp7 = ( d - c - h1/sigma*(u1+tmp0) ) * s2
         h2 = ( tmp6*tmp2/s1 - p2*tmp7 ) / d1
         h3 = - ( tmp6*tmp3*s1 - p1*tmp7 ) / d1
         h4 = -f*p3 - c*d
         tmp8 = h4/sigma
         tmp9 = ( u3 - tmp8*(u2-tmp0) ) * s2
         h5 = - ( tmp8*tmp4/s1 + tmp9 ) / d2
         h6 = ( tmp8*tmp5*s1 + tmp9 ) / d2
         h7 = (c*tmp2) / (d1*s1)
         h8 = (-c*tmp3*s1) / d1
         h9 = tmp4 / (d2*s1)
         h10 = (-tmp5*s1) / d2

!downward direct and diffuse fluxes below vegetation

         if (ic .eq. 0) then
            ftds = s2
            ftis = h4*s2/sigma + h5*s1 + h6/s1
         else
            ftds = 0.
            ftis = h9*s1 + h10/s1
         end if
         ftd(ib,i) = ftds
         fti(ib,i) = ftis

!flux reflected by vegetation
 
         if (ic .eq. 0) then
            fres = h1/sigma + h2 + h3
         else
            fres = h7 + h8
         end if
         fre(ib,i) = fres

!flux absorbed by vegetation

         fab(ib,i) = 1. - fre(ib,i) - (1.-albsod(ib))*ftd(ib,i)& 
     &               - (1.-albsoi(ib))*fti(ib,i)

      end do

      return
      end


!----------------------------------------------------------------------------------------
!apar simulations for stomatal resistance estimate by LSM_Bonan Optinal, added by AYu
!-----------------------------------------------------------------------------------------
     subroutine surrad (fsun, elai, esai, vai, solad,solai, fabd, fabi, ftdd, &
     & ftid, ftii, albsod, albsoi, laisun, laisha, parsun, parsha)
!------------------------ code history ---------------------------
!source file:       surrad.F
!purpose:           solar fluxes absorbed by vegetation and ground surface
!date last revised: March 1996 - lsm version 1
!author:            Gordon Bonan
!standardized:      J. Truesdale, Feb. 1996
!reviewed:          G. Bonan, Feb. 1996
!-----------------------------------------------------------------
!------------------------ input/output variables -----------------
!input

      use INMSOIL_NUMPARAMS, only : &
     & ML, MS

       INCLUDE 'soil.inc'
!      PARAMETER (nv2=13) 
!      USE PTRMODPHY, ONLY : NV2 

      PARAMETER (twoband=2)
      REAL (KIND=8)  fsun(nv2)          !sunlit fraction of canopy
      REAL (KIND=8)  elai(nv2)          !leaf area, one-sided
      REAL (KIND=8)  esai(nv2) 
      REAL (KIND=8)  vai(nv2)           !leaf + stem area, one-sided

      REAL (KIND=8)  solad(twoband)       !incodmin1g direct solar radiation (w/m**2) 
      REAL (KIND=8)  solai(twoband)       !incodmin1g diffuse solar radiation (w/m**2) 
      REAL (KIND=8)  fabd(twoband,nv2)    !flux abs by veg (per unit incodmin1g direct flux) 
      REAL (KIND=8)  fabi(twoband,nv2)    !flux abs by veg (per unit incodmin1g diffuse flux)
      REAL (KIND=8)  ftdd(twoband,nv2)    !down dir flux below veg (per incodmin1g dir flux) 
      REAL (KIND=8)  ftid(twoband,nv2)    !down dif flux below veg (per incodmin1g dir flux) 
      REAL (KIND=8)  ftii(twoband,nv2)    !down dif flux below veg (per incodmin1g dif flux) 
      REAL (KIND=8)  albsod(twoband)      !ground albedo (direct)
      REAL (KIND=8)  albsoi(twoband)      !ground albedo (diffuse)
!     REAL (KIND=8)  albd(twoband,nv2)    !overall surface albedo (direct)
!     REAL (KIND=8)  albi(twoband,nv2)    !overall surface albedo (diffuse)

!output      
      REAL (KIND=8)  laisun(nv2)        !sunlit leaf area index, one-sided
      REAL (KIND=8)  laisha(nv2)        !shaded leaf area index, one-sided
      REAL (KIND=8)  parsun(nv2)        !average absorbed par for sunlit leaves (w/m**2)
      REAL (KIND=8)  parsha(nv2)        !average absorbed par for shaded leaves (w/m**2)
!      REAL (KIND=8)  ndvi(nv2)          !normalized difference vegetation index
!-----------------------------------------------------------------

!------------------------ local variables ------------------------
      integer i               !loop counter/array index
      integer ib              !waveband number (1=vis, 2=nir)

      REAL (KIND=8)  abxx               !absorbed solar radiation (w/m**2) 
      REAL (KIND=8)  rnir               !reflected solar radiation [nir] (w/m**2)
      REAL (KIND=8)  rvis               !reflected solar radiation [vis] (w/m**2)
      REAL (KIND=8)  laifra             !leaf area fraction of canopy
      REAL (KIND=8)  trd                !transmitted solar radiation: direct (w/m**2)
      REAL (KIND=8)  tri                !transmitted solar radiation: diffuse (w/m**2)
      REAL (KIND=8)  cad(twoband,nv2)     !direct beam absorbed by canopy (w/m**2)
      REAL (KIND=8)  cai(twoband,nv2)     !diffuse radiation absorbed by canopy (w/m**2)
	  REAL (KIND=8)  fsha(nv2)          !shaded fraction of canopy
	  REAL (KIND=8)  sav(nv2)           !solar radiation absorbed by vegetation (w/m**2)
      REAL (KIND=8)  sag(nv2)           !solar radiation absorbed by ground (w/m**2)
      REAL (KIND=8)  fsa(nv2)           !total absorbed solar radiation (w/m**2)
      REAL (KIND=8)  fsr(nv2)           !total reflected solar radiation (w/m**2)
      REAL (KIND=8)  mpe                !prevents underflow errors if division by zero

!-----------------------------------------------------------------

	mpe = 0.000001
	do ib = 1,twoband
		do i = 1,nv2
!			if(albd(ib,i).eq.0.) albd(ib,i) = 0.2
!			if(albi(ib,i).eq.0.) albi(ib,i) = 0.2
			if(fabd(ib,i).eq.0.) fabd(ib,i) = 0.2
			if(fabi(ib,i).eq.0.) fabi(ib,i) = 0.2
			if(ftdd(ib,i).eq.0.) ftdd(ib,i) = 0.2
			if(ftid(ib,i).eq.0.) ftid(ib,i) = 0.2
			if(ftii(ib,i).eq.0.) ftii(ib,i) = 0.2
		end do
		if(albsod(ib).eq.0.) albsod(ib) = 0.2
		if(albsoi(ib).eq.0.) albsoi(ib) = 0.2
	end do
!zero summed solar fluxes

      do i = 1, nv2
         sag(i) = 0.
         sav(i) = 0.
         fsa(i) = 0.
      end do

!loop over twoband wavebands

      do ib = 1, twoband
         do i = 1, nv2

!absorbed by canopy

            cad(ib,i) = solad(ib)*fabd(ib,i)
            cai(ib,i) = solai(ib)*fabi(ib,i)
            sav(i) = sav(i) + cad(ib,i) + cai(ib,i)
            fsa(i) = fsa(i) + cad(ib,i) + cai(ib,i)

!transmitted = solar fluxes incident on ground

            trd = solad(ib)*ftdd(ib,i)
            tri = solad(ib)*ftid(ib,i) + solai(ib)*ftii(ib,i)

!solar radiation absorbed by ground surface

            abxx = trd*(1.-albsod(ib)) + tri*(1.-albsoi(ib)) 
            sag(i) = sag(i) + abxx
            fsa(i) = fsa(i) + abxx
         end do
      end do

      do i = 1, nv2

!partion visible canopy absorption to sunlit and shaded fractions
!to get average absorbed par for sunlit and shaded leaves

!	   vai(i) = elai(i) + esai(i)
         fsha(i) = 1.-fsun(i) 
         laisun(i) = elai(i)*fsun(i)
         laisha(i) = elai(i)*fsha(i)
          ! write(*,*) i, fsun(i)
         laifra = elai(i) / dmax1(vai(i),mpe)
         if (fsun(i) .gt. 0.) then
            parsun(i) = (cad(1,i)+fsun(i)*cai(1,i)) * laifra / &
     &                  dmax1(laisun(i),mpe)
            parsha(i) = (fsha(i)*cai(1,i))*laifra / dmax1(laisha(i),mpe)
         else
            parsun(i) = 0. 
            parsha(i) = (cad(1,i)+cai(1,i))*laifra /dmax1(laisha(i),mpe)
         endif

!ndvi and reflected solar radiation

!         rvis = albd(1,i)*solad(1) + albi(1,i)*solai(1) 
!         rnir = albd(2,i)*solad(2) + albi(2,i)*solai(2)
!         fsr(i) = rvis + rnir
!         ndvi(i) = (rnir-rvis) / dmax1(rnir+rvis,mpe)
      end do

      return
      end
!------

